function [Q,sphi,empty] = drawQ0_ar(restr,phi,opt)
% Function attempts to draw Q from the space of orthonormal matrices 
% satisfying sign, elasticity and narrative restrictions.
% Inputs:
% - restr: structure containing information about restrictions
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

signRestr = restr.signRestr;
elastRestr = restr.elastRestr;
shockSignRestr = restr.shockSignRestr;
hdSignRestr = restr.hdSignRestr;
Sigmatr = phi.Sigmatr;
Sigmatrinv = phi.Sigmatrinv;
vma = phi.vma;
U = Sigmatrinv*restr.U; % Multiply by Sigmatrinv to avoid doing this repeatedly
Sigmatrinvp = Sigmatrinv'; % Transpose to avoid doing this repeatedly

n = size(Sigmatrinv,1); % Number of variables in VAR
mS = size(signRestr,1); % No. of traditional sign restrictions
mE = size(elastRestr,1); % No. of elasticity restrictions
mSS = size(shockSignRestr,1); % No. of shock-sign restrictions
mHD = size(hdSignRestr,1); % No. of restrictions on historical decomp.

%% Construct matrix representing traditional sign restrictions.
% Restrictions are represented as S(phi)*vec(Q) >= 0.

S = zeros(mS,n^2);

for ii = 1:mS % For each restriction
    
    if signRestr(ii,5) == 1 % Sign restriction on impulse response
    
        S(ii,(signRestr(ii,2)-1)*n+1:signRestr(ii,2)*n) = ...
            vma(signRestr(ii,1),:,signRestr(ii,3)+1)*signRestr(ii,4);
    
    elseif signRestr(ii,5) == 2 % Sign restriction on A0
        
        S(ii,(signRestr(ii,2)-1)*n+1:signRestr(ii,2)*n) = ...
            Sigmatrinv(:,signRestr(ii,2))'*signRestr(ii,4);
        
    end
    
end

%% Construct matrix representing elasticity restrictons.
% Restrictions are represented as ES(phi,m)*vec(Q) >= 0, where m is a 
% vector containing the bounds. Do this by writing restrictions as 
% m*B(k,l) - s*B(i,j) >= 0.
ES = zeros(mE,n^2);

for ii = 1:mE

    ES(ii,(elastRestr(ii,4)-1)*n+1:elastRestr(ii,4)*n) = ...
        elastRestr(ii,6)*Sigmatr(elastRestr(ii,3),:); % m*B(k,l)
    ES(ii,(elastRestr(ii,2)-1)*n+1:elastRestr(ii,2)*n) = ES(ii,(elastRestr(ii,2)-1)*n+1:elastRestr(ii,2)*n) ...
        -elastRestr(ii,5)*Sigmatr(elastRestr(ii,1),:); % -s*B(i,j) - first term allows for possibility that j = l

end

%% Construct matrix representing shock-sign restrictions.
% Restrictions are represented as SS(phi,U)*vec(Q) >= 0.
SS = zeros(mSS,n^2);

for ii = 1:mSS % For each restriction
    
    SS(ii,(shockSignRestr(ii,1)-1)*n+1:shockSignRestr(ii,1)*n) = ...
        U(:,shockSignRestr(ii,3))'*shockSignRestr(ii,2);
    
end

% Combine coefficient matrices from restrictions.
sphi = [S; ES; SS];

%% Attempt to draw Q satisfying restrictions.

iter = 0;
empty = 1;

while iter <= opt.qMax && empty == 1
    
    % Draw Q from space of orthonormal matrices.
    z = randn(n);
    [Qtilde,~] = qr(z);
    
    % Normalise diagonal elements of A0 to be positive. Note that Matlab 
    % is implicitly expanding arrays to be compatible with elementwise 
    % array operations.
    Q = ((sign(diag(Sigmatrinvp*Qtilde))').*Qtilde)./vecnorm(Qtilde);
    
    % Check whether proposed draw satisfies sign, elasticity and shock-sign
    % restrictions.
    sSat = all(sphi*Q(:) >= 0);
       
    if sSat == 0  
        % If restrictions not satisfied, increment iteration
        % counter and return to beginning of loop.
        iter = iter + 1;
        continue;          
    end

    % Check whether proposed draw additionally satisfies restrictions on 
    % historical decomposition.
    hdCheck = zeros(mHD,1); 

    for ii = 1:mHD % For each restriction  
    
        % Pre-compute q_j*q_j' for j=1,...,n.
        QQ = pagemtimes(reshape(Q,[n,1,n]),reshape(Q,[1,n,n]));

        % Extract row of hdSignRestr corresponding to ith restriction.
        hdRestr = hdSignRestr(ii,:);

        % Compute contribution of each shock.
        HH = zeros(1,1,n);    

        for hh = 1:hdSignRestr(ii,4)+1 % Sum contribution over horizons

            HH = HH + pagemtimes(pagemtimes(vma(hdRestr(1),:,hh),QQ),...
                U(:,hdRestr(3)+hh-1));

        end

        if hdRestr(5) == 1 && hdRestr(6) == 1

            % Type A - most important contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) == max(abs(HH));

        elseif hdRestr(5) == 1 && hdRestr(6) == -1

            % Type A - least important contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) == min(abs(HH));

        elseif hdRestr(5) == 2 && hdRestr(6) == 1

            % Type B - overwhelming contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) - ...
                sum(abs(HH((1:n) ~= hdRestr(2)))) >= 0;                

        elseif hdRestr(5) == 2 && hdRestr(6) == -1

            % Type B - negligible contributor
            hdCheck(ii) = abs(HH(hdRestr(2))) - ...
                sum(abs(HH((1:n) ~= hdRestr(2)))) <= 0;

        end
        
        if hdCheck(ii) == 0              
            % If ith restriction on historical decomposition not 
            % satisfied, exit loop (saves computing contributions 
            % relating to subsequent restrictions).
            continue;
        end

    end

    hdSat = all(hdCheck);
    
    % If all restrictions satisfied, terminate while loop. Otherwise,
    % increment counter and attempt to draw again.
    empty = 1 - sSat*hdSat;
    iter = iter + 1;
        
end

end